Thermal equilibration between two quantum systems 
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Two identical finite quantum systems prepared initially at different temperatures, isolated from 
the environment, and subsequently brought into contact are demonstrated to relax towards Gibbs- 
like quasi-equilibrium states with a common temperature and small fluctuations around the time- 
averaged expectation values of generic observables. The temporal thermalization process proceeds 
via a chain of intermediate Gibbs-like states. We specify the conditions under which this scenario 
occurs and corroborate the quantum equilibration with two different models. 
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The derivation of thermodynamic phenomena from de- 
terministic time-reversible dynamics constitutes one of 
the primary goals of physics. This long-standing conun- 
drum has sparked recently a new wave of activity in the 
quantum domain, where current studies of the objective 
follow essentially two tracks. The first one, pioneered 
by Schrodinger [1], leads to an understanding of canoni- 
cal thermalization when the system of interest is coupled 
to a much larger system, a quantum Giant [2-6]. The 
studies along the second track explore the "microcanoni- 
cal" thermalization within a single isolated quantum sys- 
tem [?]-[l2], and employ exact numerical diagonalization 
of many-body models [13l [l4 1 . 

Here, we focus on a different route by elucidating the 
process of mutual equilibration between two finite quan- 
tum "peers" , prepared initially at different temperatures 
and then set into a contact. We consider two systems, 
A and B, that are identical, in the sense that they have 
identical Hamiltonians, Ha = Hb = Hg. The Hamil- 
tonian Hs has J\fs non-degenerate energy levels, {e/c}, 
k = l,...,J\fs, i.e., H s \4*k) = e fe |0fe), with eigenstates 
{!</>*;)}■ The systems interact through a contact, which 
allows only for energy transfer without exchange of parti- 
cles. The Hamiltonian of the composite bipartite system 
thus reads 



H< 



XH 



int. 



(1) 



with A being a dimensionless coupling constant. The 
interaction Hamiltonian, H mt = Ya ® Yb, with operators 
Ya = Yb = Y, is invariant under permutation 4 f> B 
and does not commute with the Hamiltonian Hs [lij . 

We denote the energy eigenvalues and the corre- 
sponding eigenstates of the Hamiltonian H x by {E x } 
and respectively. The quantities of interest, 

i.e. the energy level populations, p A (t) and pf(t), can 
conveniently be calculated by using the product basis, 
l^n(fe j)) = ® which is also the eigenbasis of the 
composite system for A = 0. We label the energies E® ac- 
cording to their decomposition into the sum of the single 
system energies, E®^ = f-k + £j = E n (j,k)- To shorten 
notations, we shall use either n or kj instead of n(k,j). 
While combinations k — j produce the non-degenerate 



energy levels, E kk = 2ek, each two levels related by the 
permutation of indices k •<-> j, with k ^ j, are doubly de- 
generate; i.e. E^j = Ej k . The transformation from the 
product basis \ip%) to the eigenbasis at a certain interac- 
tion strength A > 0, \ipn), is given by the matrix A, with 
the elements A„ im = (ipm\4>n) ■ Throughout this work we 
further assume for the Hamiltonian ([l} with A ^ both 
the non-degeneracy, E* ^ E^ for n ^ m, and the "non- 
degenerate energy gap condition" 0, [f| [HI , meaning 
that non-zero energy differences E* — E^ and E* — E x 
are not equal, apart from the trivial case s — n, w = m. 

The energy level populations p A (t) for system A{B) 
are given by the partial trace over system B(A) of 
the composite system density matrix g(t); for example, 
Pk M = J2j Qkj,kj(t), where g(t) is expressed in the prod- 
uct basis. In the case of canonical initial states, where 
only diagonal density matrix elements are initially non- 
zero, their evolution can be described by the linear map 



Qn,n{t) 



(2) 



where U*(t) 
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It is apparent 



that all necessary information is encoded in the energy 
spectrum {E x } and in the transformation matrix A. 

For any choice of the system initial states, g A (0) and 
g B (0) , the mutual equilibration is guaranteed (in a sense 
detailed below) as long as the non-degenerate energy gap 
condition holds. Due to the parity A ■<-» B all eigen- 
states of the Hamiltonian H x are either symmetric, \tp k j) 

— or antisymmetric IV'jy) = ~\' l Pjk)- Therefore, 

for every eigenstate of the composite system, expecta- 
tion values for any local observable O (energy, level pop- 
ulations, etc.), associated with one quantum peer only, 
would be the same for the second peer, A = B . Hav- 
ing the total system prepared at time t = in a prod- 
uct state g(0) — g A (0) ® g B {0), we turn on the inter- 
action by setting A > 0. Then, after some characteris- 
tic relaxation time r re /, the system is expected to reach 
quasi-equilibrium, where all diagonal elements of the two 
subsystem reduced density matrices obey the relation 
g kk (t) ~ Qkk{t) 0- The respective total equilibrium 
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system energies can be evaluated from the condition of 
energy conservation (assuming a diminutive interaction 
energy), 



E. 



A.B 



\[E A (0) 



E B (0)], 



(3) 



where E° — J2k e kQk k w hh S = A or B. This is not a 
genuine equilibrium, since the populations still evolve in 
time (l7| . but their recurrences occur on time scale T rec 
which is larger than any relevant time scale 12|, [18J . 

To gain an analytical insight, we start out from the 
limiting case in which the transformation matrix takes 
on a simple form: Any infinitesimally small interaction, 
A — s- 0, will lift the two- fold degeneracy, E^ = E^ k , 
yielding the pair of a symmetric and an antisymmetric 

eigenstates in the form (l^kj) ^ l^jtS) > ^ 3- These 
eigenstates are non-degenerate and separated by a finite 
splitting. The eigenstates whose energies, E® k = 2e k , 
were non-degenerate at A = 0, are perturbed marginally 
only in this limit. By assuming this so resulting tridiago- 
nal structure for the transformation matrix A„ jTO , we find 
that the relaxation process leads to the arithmetic-mean 
quasi- equilibrium state, with the corresponding popula- 
tions reading [19( 



A.B 

Pk 



[P*(0)+Pf(0)] 



(4) 



This tridiagonal structure is guaranteed to hold as long as 
each off-diagonal, non-zero matrix element of the inter- 
action Hamiltonian, AI-H^J, is smaller than the corre- 
sponding energy level difference in the composite system, 
AE n , m = \E° n -E° m \. 

The characteristic feature of the arithmetic- mean equi- 
libration is that two systems, when initially prepared in 
canonical states at different temperatures, gf an (Ts) — 
e -H s /k B T s j Zs ^ z s = Tr(e- Hs / kBTs ), with the diagonal 
elements 



glk(T s )=P S k 



1 p -e k /k B Ts 



(5) 



where fcg is the Boltzmann constant, do relax to states 
with the same mean energy, but their energy level pop- 
ulations, Eq. ((3]), are no longer Gibbs-like. In order to 
deviate from the limit in Eq. Q the transformation ma- 
trix A needs to acquire a more complex structure. This 
is achieved by cranking up the interaction strength be- 
tween the two systems. Provided that there occurs a 
sufficiently large number of non- vanishing off-diagonal el- 
ements, increasing the strength of interaction, but 
still remaining within the weak coupling limit 



_int 



) < ejv s 
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wherein {e™*} is the spectrum of the interaction Hamil- 
tonian H lnt , then yields interaction blocks in the ma- 
trix A n?m larger than those 2x2 blocks. We expect 
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FIG. 1: (color online) (a) A system of bosons confined into two 
overlapping confinements is analyzed with the Bose-Hubbrad 
model, (b) Energy spectrum of a single system. The (red) 
line displays the dependence of the system mean energy, i.e., 
E s = J2k £ke~ ek/kBT /Z s , on temperature T. The initial tem- 
peratures of the 'hot' system, fceT^/s = 94.91, and the 'cold' 
system, IcbTb/s — 18.98, are indicated by the (blue) dots. 
The equilibrium temperature, feeTp/s = 33.92, calculated by 
using the total energy conservation, Eq. (J7J), is indicated by 
the (red) star, (c) Instantaneous 'equilibrium' energy level 
populations for systems A (left column) and B (right column) , 
in the regime of arithmetic-mean (top) and thermal (bottom) 
equilibrations. The arithmetic-mean populations are depicted 
by the top (blue) solid lines, and the canonical populations for 
the temperature Tf by the bottom (red) lines. The natural 
energy unit, s, is given by the mean energy level spacing of the 
single system, s = (ejv g ~ £i)/(A/"s — 1). The similar behav- 
ior is demonstrated by the second model, see supplementary 
material for further model details. 



that the presence of a more complex block structure en- 
sures the evolution of canonical initial states, g£ an (TA) 
and gf an (Tb ) , towards a common Gibbs-like equilibrium 
q a ' b (Tf), meaning that the corresponding diagonal el- 
ements are given by the relation ([5]) with the common 
temperature Tp. The 'equilibrium' temperature Tp can 
be evaluated from Eq. ([3]), to yield with Eq. 



Z F 



e k B T A 
Z A 
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(7) 



We numerically validate our prediction by using two 
types of quantum models. Within the Bose-Hubbard 
model we consider the system consisting of N = 5 on- 
site interacting bosons on a one-dimensional lattice, with 
L = 5 sites and hard-wall boundaries. This results in 
■^s = fe~TTTOT — 126 energy levels in each single sys- 
tem, and JV — Ms x As = 15, 876 levels in the com- 
posite system [ijj]. Figure 1(a) depicts the setup, which 
assumes that the two systems overlap only by one site, 
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FIG. 2: (color online) Relaxation pathways for the model 
depicted in Fig. 1(a). Both systems are initially prepared 
in canonical states (solid lines) and in pure states randomly 
sampled from the corresponding ensembles of typical states, 
Eq. ((SJ (dashed lines), (a) The evolution of the mean energies 
and (b) the corresponding temperatures T(t) of the 'hot' sys- 
tem A and the 'cold' system B are shown by the upper (red) 
and lower (blue) lines, respectively, (c) The energy level pop- 
ulations of both systems are displayed at different moments of 
time (dots), marked by the corresponding symbols in (a, b). 
The lines correspond to the canonical populations, Eq. ([5J, at 
the temperatures evaluated from the temporal values of mean 
system energies (see Fig. 1(b)). 



where the bosons from the different confinements do in- 
teract. We also corroborated our findings with a ran- 
domly synthesized model, for which the Hamiltonian H$ 
and the interaction operator Y are independently sam- 
pled from a finite-dimensional Gaussian Orthogonal En- 
semble (GOE) of random matrices [13 ], In contrast to 
the former many-body interacting boson model, where 
the interaction is strictly local, here the interaction is 
now acting globally, interweaving systems A and B. 

For both models we find solutions that are based on 
the exact diagonalization of the corresponding bipartite 
Hamiltonians. Our main results are depicted in Fig. [1] 
Upon increasing the coupling constant A within the weak 
coupling limit, Eq. ©, we detect a crossover from the 
arithmetic- mean quasi-equilibrium populations, Eq. ((][]), 
towards the canonical populations, Eq. ([5]) with T5 = Tp. 

An intriguing question is how the quantum equilibra- 
tion unfolds in time. Figure 2 displays our finding that 
equilibration proceeds along a quasistatic pathway: the 
relaxation of an initial canonical state abides a sequence 
of time-dependent Gibbs-like states with time-dependent 
temperatures T(i), intermediate between the initial tem- 
perature T4( B ), to reach a common, final temperature 
Tp. This observed persistence of Gibbs shape is remark- 
able indeed. The only relevant result we could find in 
this context is that of thermal relaxation dynamics of a 
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FIG. 3: (color online) (a) Von Neumann entropy of a single 
system vs. time for the model shown in Fig. 1(a). Both sys- 
tems are initially prepared in pure states randomly sampled 
from the ensembles of typical states, Eq. JHJ. The dashed 
line indicates the entropy of the Gibbs state at the tempera- 
ture Tf- (b) The population dynamics n?s(t) of the fifth-site 
(i.e. the site making the thermal contact) for the subsys- 
tem A is compared with the corresponding canonical value 
at the temperature Tf. Note that the time average of n$(t), 
0.9673, differs from its canonical value, 0.9651, by 0.3% only, 
(c) The absolute values of reduced density matrix elements 
{4>m\Q A (t)\(f>„) at t = 0, and (d) after equilibration. 



stylized model [20[. 

We next consider the case with an initial preparation 
given by pure states. Reproducibility of quantum ther- 
mal processes with a single 'typical' state [21| carries im- 
portance in view of the foundations of statistical physics 
7, 22] and many-body quantum calculations [23J. We 
employ here typical states constructed as the sums over 
eigenstates [2l|; i.e., we use 



IV&(0)> 



1 



k 



e k /2k B T s 



I'M- 



(8) 



The ensemble of typical states is defined by the uni- 
form measure on the torus Of (g) #f ... <g> Ofr, 9§ £ [0,27r]. 
The results shown in Fig. 2 (a, b) by the dashed lines 
confirm our expectation: A single, randomly sampled, 
initial product wave function |^(0)) <g> |V# B (0)) fol- 
lows the equilibration pathway for canonical initial states 
with good accuracy. Both systems, A and B, are pre- 
pared initially in pure states, implying vanishing von 
Neumann entropies S A ,B(t) = -k B Tr[g A ^ B (t) In g A ' B (t)}, 
i. e., Sa(0) — Sb(0) — 0. The isolated composite system 
remains in a pure state forever, and thus SA®B(t) = 0. 
This, however, is no longer so for subsystem entropies 
SU(£) an d Ssii), which start to grow. From the tri- 
angle inequality it follows that SU(i) = = S(t). 
The entropy S(t) is a measure for entanglement between 
the subsystems [2^|: its monotonic growth thus indi- 
cates that the equilibration process entangles the quan- 
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turn peers, see Fig. 3(a). 

The systems cannot rigorously reach canonical equilib- 
rium; therefore, the entropy S(t) saturates to the value 
below the entropy of the Gibbs state at temperature Tp. 
The resulting equilibrium system density matrices, q (t) 
and g B (t), remain nonstationary and possess both diago- 
nal and off-diagonal elements evolving in time. Following 
the recipe from Ref. 0], the deviation from the canoni- 
cal state is estimated by using the trace-norm distance 
V = Tr (\g s (t) — g s \)t/2, where the bar denotes the time 
average (. . . ) t . This quantity is limited from above @, 
so that from Eq. (8) in Q we find that V < 0.6 in our 
case. From our numerics we obtain T> ~ 0.43. 

For an operator O, which is non-diagonal in the eigen- 
basis of the Hamiltonian Hs, the presence of the off- 
diagonal elements in the system density matrices will 
produce additional fluctuations around the average value 
s = Tr(g s O). Moreover, some of the off-diagonal ele- 
ments may possess non-zero time averages. This might 
cause a constant shift of the observable averaged value 
O from its canonical value, SO s = s - Tr[pf an (T F )0]. 
However for highly non-sparse patterns of non-zero off- 
diagonal elements 0^ and g^ , we may expect that the re- 
spective fluctuations of the expectation value s (t) will 
be suppressed, exhibiting dynamical typicality [25j |. Even 
for a system as small as ours, with Afg — 126 states, this 
mechanism works surprisingly well, see Fig. 3(b). 

Thermal quantum relaxation within an isolated com- 
posite quantum system is a deterministic process, and 
produces an output in the form of a Gibbs-like equilib- 
rium, with diagonal elements which are almost canonical, 
for the initial preparation, Eq. ([5]), and also for initial 
'typical' pure states, Eq. ((5J). An arbitrary choice of the 
initial state of the composite system H x does not guar- 
antee relaxation towards Gibbs-like quasi-equilibrium 
states for its halves. Also the state of the composite 
system after relaxation is far from being Gibbs-like due 
to strong entanglement between its halves. Moreover, in 
order to render the thermodynamical relaxation of quan- 
tum peers, two necessary conditions need to be fulfilled, 
namely, (i) the interaction is restricted to the validity 
range of Eq. (j6|), and (ii) the total composite system 
obeys the parity A o B. A natural question then is: 
What will happen if either of the conditions (i) or (ii) is 
violated? For (i) the systems will nevertheless equilibrate 
even with the interaction strength set beyond the weak- 
coupling limit. The corresponding 'equilibrium' state, 
however, no longer assumes a Gibbs-like structure. The 
part (ii) with non-identical systems A and B is more in- 
tricate. Although it is still possible to obtain thermal 
relaxation between two different systems (see [l9]), the 
mismatch of system spectra and their relatively small 
sizes necessitates a much larger system-system coupling 
constant A . The resolution of this problem demands 
systems of much larger sizes, and, therefore, lies outside 
the exact diagonalization scheme employed here. 



The quasistatic character of the thermal relaxation al- 
lows for the tuning of one of the two quantum peers to a 
Gibbs-like state at any temperature between initial tem- 
perature values, Ta and Tg, thus serving as an alter- 
native protocol for the preparation of thermal states of 
quantum systems [26[. The state-of-the-art experiments 
with ultracold atoms provide the natural playground for 
exploration of the thermal relaxation between two differ- 
ent species of atoms [53] . 

This work is supported by the DFG grant HA1517/31- 
2 and by the "Nanosystems Initiative Munich" (NIM). 
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ARITHMETIC-MEAN EQUILIBRATION 

The arithmetic-mean equilibration occurs in the limit 
of weak coupling between the systems, \{H mt } <C {AE}, 
where the transformation matrix A n . m = (V'ml^/n) effec- 
tively takes on the tridiagonal form, see Fig. HJa), with 
non-zero entries: 

f 1, j=f = k= k'; 
K{k,j),m{k',j') = \ 1/V2, j=j'^k = k f \ (1) 
[ Xk 3 /V2, j = k'^k=f, 

where Xkj = sgn(fc — j). Hence the dynamics of the com- 
posite system is governed by the tridiagonal evolution 
operator: 



1,, 



le^Xk'j'Sin [ 



m,n^rn,ri 

J=.f = k = k>; 
j=f^k = k>; 



(2) 



, j^k,j = k'^k = f, 



where fi* = (E* 



E} k )/2h and w * = (E* 



For the density matrix of the composite system written 
in the product basis, knowledge of its diagonal 

elements is sufficient for calculations of the energy level 
populations. For example, for the system A one has 



Pk (*) = ^2 Sn(k,j),n(k,j)(t)- 



(3) 



The energy level populations of one system, e.g. the sys- 
tem A are 

Pt(t) = ~ [^(0) + Pf (0)] + W X* cos ( w *f) , (5) 



where X§ = p£(0)pf (0) - pf(0)p% (0). 

In order to gain some physical insight it is useful to ex- 
plore the following situations: (i) one of the two systems, 
for an example system B, is initially localized on a single 
level, i.e. p k = 8k,k B ^ an d (ii) both the systems, A and 
B, are initially localized on single levels, i.e. p k = 5k,k A 
and pf = S k ,k B ■ 

In the first case (i), Eq. JS]) reduces to 



E^p/(0)cosK a /), k = k B ; 
Pfe(0)cos(w A t), k + k B , 



while for the energy level populations of the system B we 
obtain, 



Pk(t) 



+ Pfe(0) + £, ¥fc p/(0)cos(^), 

+p£(o) 



Pfc(0)cos« s i), M*i 



Thus the population of the energy levels with the same 
index k 7^ k B in the system A and the system B oscillate 
coherently, while the populations of energy levels with 
k = k B quasiperiodically fluctuate, for all k around the 
arithmetic-mean of the initial populations. 

In the second case (ii), the energy level populations of 
systems A and B read: 



Evolution of initially diagonal states 



For the systems whose initial density matrices are di- 
agonal in {|0°}}, p A ' B (0)k,k' = 4,fc'Pfc" B (0), the evolu- 
tion of the diagonal matrix elements for the composite 
system density matrix, p n (t) = g n ,n(t), reduces to the 

linear map: p n (t) = J2 n > l^n.n'Wl Pn'ity- % usin S the 
evolution operator, Eq. (J2|), we obtain 



Pkj(t) 



(p kj (0)+p jk (0)) 



(^)(p fej (0)-p ife (0)). (4) 



Pkit) 



and 




Pk(t) 



cos w 




( w fcfc H 0, k = k A ; 



C0S ( w Lfc*)i k = k B ; 



k ^ k A ,k ^ k B ; 



cos ( w fcfc f 

C0S ( w fe A fe*)> k = k B 



t)i k — k A ; 
k = k B ; 

k ^ k A ,k ^ k B . 



Thus, for k A 7^ k B each system performs coherent Rabi 
type oscillations between the levels k A and k B around 
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FIG. 1: Structure of transformation matrix, (a) The part of transformation matrix, A n n /, for the GOE model of quantum 
peers in the limit of arithmetic-mean equilibration, (b) The same part of transformation matrix beyond the arithmetic-mean 
limit, (c) The model of the 36 x 36 transformation matrix with four 2x2 blocks, four 4x4 blocks, one 6x6 block, and the 
seven corresponding free parameters, at, i = 1, 7. Only the absolute values of transformation matrix elements are shown. 



the arithmetic-mean of their initial values, while the rest 
of the levels remains unpopulated. Evidently, if k^ — 
kg = kg then the corresponding frequency uJk s k s = 0> 
thus both systems stay localized forever and there is no 
dynamics at all. 



Evolution of initially pure states 



When initially both the systems are assumed to be 
in pure states, \tp A,B (0)) = Xlfc c fc '^l^fc) an< ^ the initial 
state of the composite system is given by their product, 
1^(0)) = \ip A (0)) ® \ip B (0)), so that 



W°)> = 5Z c «(fe,i)IV'°(fc J )>> c n = c A c 



A„B 



(6) 



the diagonal matrix elements of the total density matrix 
are 



p kj (t) = | «#(*)> | a = |(V^-|u A (t)|^(0))| 2 . 

Of)} , where 9 A 



(7) 



i /> 



Using c k j = \JP A pf exp [j 
the phase entering the initial state of the system A (B) 
(see Eq. (8), of main manuscript), we end up with the 
the energy level populations of the system A: 

p A (t) - \ [p A {Q)+p B {Q)] + l -Y. x h co « t ) 

3 

+\ E y 4 sin - sin «*) » ( g ) 



where Y A . 
6 A , as 



= X J fcV / Pfc(0)pj J (0)p/(0)pf(0)/2 and 6 kj = 

,B . The only difference from the Eq. ([5]) is the 
last sum on the rhs. The latter is generated by non-zero 
off-diagonal elements of initial density matrix, g(0) = 



Apparently, the evolution for diagonal initial states can 
also be obtained from Eq. (J5J by averaging over the en- 
semble of pure states with the different initial phases, O^j , 
and fixed initial populations of the energy levels, p A ' B '(0). 
This would lead to the nullification of the last sum on the 
rhs of Eq. ©. 



Relaxation to the arithmetic-mean of the initial 
populations 

In most physical situations, the initial energy level pop- 
ulations, p A ' B (0), exhibit a smooth dependence on k, 
thus producing a significant number of non-zero coeffi- 
cients X A j B . This is also the case for the canonical states 
with fcgT 3> s, where s denotes the mean level spacing, 
s = ( e .V s - ei)/(A/"s - 1). The relation (J5J yields p A (0) 
at t = 0. In the course of time every member of the sum 
on the rhs of Eq. ([5]) begins acquiring a certain phase. 
Since the frequencies do not commensurate in gen- 
eral, after the characteristic time, T re \ 



2ir/ujt yp , where 

uu typ is the root mean square of the set {wjy}, we will 
obtain a sum of independent random values, almost uni- 
formly distributed over the interval [—1, 1], and weighted 
with the coefficients X A y Given the initial states with a 
substantial number of populated energy levels, the sum 
looses its initial coherence completely after the time r re ; 
and averages itself to zero. Hence the relaxation process 
leads to the arithmetic-mean equilibration, 



A,B 

Pk 



p A (0)+p B (0)] 



(9) 



The rhs of Eq. ([5]) is a quasiperiodic function Q , there- 
fore it will repeat itself after some time r rec with any 
given accuracy A, so that ||f>j^(f +T(A)) —p A (t)\\ < A. 



Assuming, for example, that the coefficients X A - B are 
equal, the recurrence time grows exponentially with Ms'- 



3 



Tree ~ -rnbsrTWs + iXAMsM-M-W i. Al- 

ready for Afg — 10 and A = 0.1 we find a sharp scale 
separation between the recurrence and relaxation times, 

^Vec/^Ve/ — 3 • 10 . 

BLOCKS IN THE TRANSFORMATION MATRIX 
A AS INITIATORS OF THERMAL 
EQUILIBRATION 

Let us denote by \tpsa) the pah, the symmetric and the 
antisymmetric eigenstates, and by ^) their two-fold 
degenerate parental eigenstates, correspondingly. Here 
index A stands for the product state, |V'Si) = |<^fe) ® |</>z), 
whereof the larger part of the energy, = e A +ef , is lo- 
cated in the system A, i.e. e A — e&, e B = e;, e& > e;. The 
eigenstate \tfj B ) is given by the permutation, \<f)i)®\(j)k) ■ In 
the limit of the arithmetic-mean equilibration the tridi- 
agonal structure of the transformation matrix allows for 
energy exchange between levels with the same energies, 
e A = e B , only. A strengthening of the interaction violates 
this 'level-to-level' rule, so that more than two eigenstates 
of the composite system can exchange their energies. Yet 
the interaction of any strength preserves the permutation 
symmetry. Therefore, a step beyond the arithmetic-mean 
equilibration consists in inclusion of a 'coupling' between 
the pair of states |Vv ) an d its neighbor (closest in to- 
tal energy) pair, \tp^ a ). In terms of the transformation 
matrix elements the coupling means that the pair g) 
contributes to the states \ip^ a ), and the pair B ) con- 
tributes to the states \ipg a )- The block corresponding to 
this unitary transformation has the following form 
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-hj 




J 



(10) 



Assuming that only one 4x4 block is present in the trans- 
formation matrix, after tracing over one system and ne- 
glecting the oscillating elements we obtain the following 
'equilibrium' populations of the energy levels: 



p£ 


= vl = 


1 

2 


[p A (0)+p B (0)]- 


SP, 


P f 


=p? = 


1 

2 


[^(o)+pf(o)]- 


SP, 


pf 


=pf = 


1 

2 


W(0)+lf(0)] + 


SP, 


pf 


=pf = 


1 

2 


[pf(0)+pf(0)]4 


-SP, 



where SP = a{p A (0)p? (0) + pf(0)pf (Q) -p6(0)pf (0) - 
pf{0)p B (0)). The only parameter is a — a 4 + c 4 + e 4 + 
.g 4 - 1/2 = b 4 + d 4 + / 4 + h 4 - 1/2, -1/2 < a < 0. 
Therefore, the 4x4 block produces the unidirectional 
energy exchange between the pairs |Vv „) and |^>* a ). 



The above result can be generalized to the case of a 
2M x 2M block, which involves the energy exchange be- 
tween M pairs of eigenstates. The energy exchange is 
then parameterized by K parameters a q , q = 1,..,K, 
where K is given by the number of ways to choose a pair 
from M elements, i.e. K = M(M - l)/2. 

When the transformation matrix has a multi-block 
structure, a given single system eigenstate \<fik) may enter 
several blocks. Then the quasi-equilibrium energy level 
population for the state \<pk} is given by 

pt B * \ [pi(o)+Pk(o)] + ( 12 ) 

where the sum index, s, runs over the numbers of those 
blocks {s}k that the fc-th state participates in, K s = 
M S (M S — 1)/ 2 is the number of parameters a q in the s-th 
block of size 2M S x 2M S , l s is the index of a partner state 
for the s-th block, and k q l q denotes the product state 
which exchanges the energy with the product state kl s . 
The typical multi-block structure of the transformation 
matrix A is depicted in Fig.QJb). 

To demonstrate that even a small number of 2M x 2M 
blocks with M > 1 in the transformation matrix A, en- 
sures that the energy level populations, Eq. ([12")) , ap- 
proach the Gibbs-like equilibrium, we use a simple model 
system with As = 6 energy levels. We take the transfor- 
mation matrix with the structure displayed in Fig. [He). 
According to (|12|). the latter is fully described by the 
seven parameters: a\ (1st 2x2 block), a\ (2nd 2x2 
block), a\,a\,a\ (3x3 block), a\ (4th 2x2 block) and 
a\ (5th 2x2 block). In the following, we rename, for sim- 
plicity, these seven parameters again to Ui, i = 1, . . . , 7. 
The Af x Af matrix, with Af = Afs x Afs = 36 entries, 
describes the composite system. The seven correspond- 
ing 'exchange' terms, <E>, composed according to Eq. (|12p . 
have the following forms: 

<*i, $u, 2 3 = p A (0)p B (0)+p A (0)p?(0)~p A (0)p!(0) 
-pf(0)p B (0), 

aa, $15,24 = pf(0)pi(0)+p A (0)p B (0)-p A (0)p B (0) 
-pf(0)pi(0), 

OS, $16,25 = P f(0)pi(0)+p A (0)pf(0)-p A (0)pi(0) 

-P5>K(o), 

«4, $16,34 = pt(0)p B (0)+pi(0)p?(Q)-pf(0)p B (Q) 
-pf(0)p 3 B (0), 

as, $25,34 = p A (0)pi(0)+pi(0)pi(0)-pf(0)p B (0) 
-pi(0)p B (0), 

«e, $26,35 = p A (Q)pi(0)+p£(0)pi(0)-p A (Q)p B (0) 
-P A (0)P$(0), 

<* 7 , $36,45 = p A (0)p B (0)+p A (0)p$(0)-pi(0)p B (0) 
-P A (0) P B (0). 
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FIG. 2: Towards the canonical populations of energy 
levels. The color diagram displays the root mean square de- 
viation, Eq. (|13[) . of the quasi-equilibrium energy level popu- 
lations from the canonical populations, for the 6-level model 
with the transformation matrix sketched in Fig. QJc). The 
three subplots (•, ♦, show the energy level populations 
(dots) in linear-log scale, for the three different sets of val- 
ues of block parameters, ot\ = otj, j = 2, . . . , 7, indicated by 
the corresponding symbols on the color diagram. In addition, 
the canonical populations (straight red lines) is plotted to- 
gether with the arithmetic-mean of initial populations (blue 
lines) . 

so the energy level populations, p^' B (0), k = 1, . . . , 6, ap- 
proach and fluctuate around the equilibrium populations 

A,B 

Pk ■ 

A B 

Pl ' = Pl + "1*14,23 + "2*15,24 + "3*16,25 + "4*16,34, 
A B 

P2 ' = PT - "1*14,23 - "2*15,24 ~ "3*16,25 

+ "5*25,34 + "6*26,35, 

A B 

P3 ' =P% "1*14,23 - "4*16,34 ~ "5*25,34 

-"6*26,35 + "7*36,45, 

A B 

PA ' = PT + "1*14,23 - "2*15,24 ~ "4*16,34 

-"5*25,34 + "7*36,45, 

A B 

P5 ' = PT + "2*15,24 - "3*16,25 + "5*25,34 

-"6*26,35 - "7*36,45, 

A B 

P6 ' = PT + "3*16,25 + "4*16,34 + "6*26,35 + "7*36,45, 

where p£ r = (p£ + p B )/2 denote the arithmetic- mean. 

Figure |2] displays the above equilibrium populations of 
the six energy levels, ej., k = 1, • . • , 6, possessed by the 
system of two bosons in three-site lattice, with the ini- 
tial temperatures fegT^ = (e2 — ei)/2, kgTs = £6 — ( i- 
To measure the deviation from the canonical population, 
we calculate the root mean square deviation (RMSD) 
of the populations pf from the canonical populations 
at the temperature Tp, evaluated from the condition 
of energy conservation, Eq. (7) of the main manuscript: 

A(ai,a 2 , • • • ," 7 ) = \/Efc=i [Pk -Pfc(T F )] 2 /6. We av- 



erage the RMSD A(ai, 07) over the region V with pos- 

A B 

itive resulting populations of the energy levels, p k ' > 0, 
in the parameter space, —0.5 < ay < 0, j = 3, . . . , 7, and, 
in Fig. [2] plot the resulting two- variable function, 

A(ai,a2)= / da^dai . . . da7A(ai, a.2, ■ ■ ■ , 0:7). (13) 
Jv 

At the point a\ = ay = 0, with j = 2, . . . , 7 (indicated by 
the symbol •, in the color diagram) the transformation 
matrix is tridiagonal and the energy level populations 
(dots) follow the arithmetic-mean of initial populations 
(solid blue lines), while any deviation from this point 
(e.g., at the two points a\ = otj, with a\ indicated by 
the symbols ♦ and -fc, in the color diagram) drags the 
energy level populations towards the canonical popula- 
tions (straight red lines, in linear- log scale). 

THE MODELS OF TWO FINITE QUANTUM 
SYSTEMS IN CONTACT 

The boson model is represented by the Bose-Hubbard 
Hamiltonian, 

J ^ U ^ 

H s = -3 X! ( a \ a i+i + a \+i a i)+^^Z n i{ n i~ l ), ( 14 ) 
1=1 1=1 

where a) (a/) is the bosonic creation (annihilation) oper- 
ator on site I, ni = ajai is the particle number operator, 
and L is the total number of lattice sites. The parame- 
ters J and U are the hopping and the on-site interaction 
strengths, respectively. The contact interaction between 
the bosons from different systems takes place at the one 
site only, Ia = L (Ib = 1), 

H int = U int n£ ®nf , (15) 

with the on-site interaction strength U mt . In all our 
calculations we set the ratio J/U = 7/3, which is far 
from the case of degenerate spectrum at J = 0, U 7^ 
or U = 0, J 5^ 0. Finally, without loss of generality, 
we set [Tint = J + U. The Hilbert space of the sys- 
tem with N bosons and L lattice sites is spanned by 
Ms = (£ + N - 1)1 /(L - l)!/iV! basis states. 

The parameters used in calculations: N = L = 5, 
J = 13.29 S, U = 5.69 s, XU int = 0.19 s (the arithmetic- 
mean equilibration, Fig. 1(c) of main manuscript, top) 
and XU lnt = 1.9 s (the thermal equilibration, Fig. 1(c) of 
main manuscript, bottom), fegT^ = 5(U + J) — 94.91 s, 
k B T B = U + J = 18.98 s, k B T F = 33.92 s. 

To synthesize the Hamiltonian of the second model, 
we employed Afg x A/j? matrix specimens from a Gaus- 
sian Orthogonal Ensemble (GOE) [l|. We use a system 
with As = 192 states, so that the composite system ex- 
hibits Af = 192 2 = 36, 864 energy levels. The matrix ele- 
ments of Hamiltonian Hs were taken from the symmetric 
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FIG. 3: Results for (a) a randomly synthesized GOE 
model of two quantum peers, (b) Energy spectrum of 
a single system. The (red) line displays the dependence of 
the system mean energy, i.e., E s = ^2 k £ke~ Ek ^ kBT /Zs, on 
temperature T. The initial temperatures of the 'hot' system, 
UbTa = 1363.81 s, and the 'cold' system, /cbTb = 34.09 s, are 
indicated by the (blue) dots. The equilibrium temperature, 
^bTf = 81.05 s, calculated by using the total energy conserva- 
tion, see Eq. (7) of main manuscript, is indicated by the (red) 
star, (c) Instantaneous 'equilibrium' energy level populations 
for the system A (left column) and B (right column), in the 
regime of arithmetic-mean (top) and thermal (bottom) equi- 
libration. The exact arithmetic-mean of initial populations is 
depicted by the top (blue) solid lines, and the canonical pop- 
ulations for the temperature Tf by the bottom (red) lines. 
The natural energy unit, s, is given by the mean energy level 
spacing of the single system, s = (e/v s ~ £i)/(Ms — 1)- 

normal distribution with dispersion a, set to one in di- 
mensionless units. It has been rescaled then to the mean 
level spacing s, so that a = 3.55s. The constraint on the 
Hamiltonian, [Hs] n ,n> = [Hs]n',n, n,n' = 1, JV, pro- 
vides its hermicity. The interaction Hamiltonian, H mi , 
was modelled by the direct product of two identical ran- 
dom matrices from GOE, generated independently, by 
using the same procedure. The parameters used in cal- 
culations are: XU lnt = 2.72- 10 -4 s (the arithmetic-mean 
equilibration, Fig. [Sfc), top) and XU int = 5.45 • 1(T 3 s 
(the thermal equilibration, Fig. [He), bottom), HbTa — 
1363.81 s, k B T B = 34.09 s, k B T F = 81.05 s. 



CROSSOVER FROM THE ARITHMETIC-MEAN 
TO THE THERMAL EQUILIBRATION 

In this section we provide the details of the crossover 
between the regimes of 'arithmetic-mean' and thermal 
equilibrations. We analyze the state of the subsystem 
A represented by its energy level populations pf{t), af- 
ter the relaxation from initial Gibbs state. To quantify 
the deviation of the actual 'equilibrium' state from the 



FIG. 4: Crossover from the arithmetic-mean to the 
thermal equilibration. The distance, Eq. {TB|, between 
the time averaged energy level populations of system A af- 
ter relaxation, and the arithmetic-mean (triangulares), and 
the thermal populations corresponding to the equilibrium en- 
ergy (dots). The crossover is shown for the three different 
combinations of particle numbers, N, and lattice sizes, L. To 
guide the eye, identical markers are connected by dashed lines, 
whose crossing points are indicated by the open circles. The 
inset shows the dependence of the single system mean level 
spacing, s = (eji/ s — e\)/{Ms — 1), Eq. (fl4)l . on the size of the 
Hilbert space, Ms, for the systems with (N, L) equal to (3, 4), 
(4,4), (4,5), (5,5), (5,6), (6,6), and (6,7). The data points 
(dots) are fitted by a power-law, s oc M$ (dashed line), with 
the exponent £ ~ —0.8. 

arithmetic-mean or Gibbs-like states, we employ the dis- 
tance between energy level populations: 

v; r >™ = Y,\pt r ' can - (p£h\ (is) 

k 

where (...)* denotes the time- average evaluated after the 
equilibration, r rec > t > T re i, and p^ r {pl an ) are the pop- 
ulations of the corresponding arithmetic-mean (Gibbs- 
like) state. 

Figure [4] shows the results of the calculations for the 
bosonic model, Eq. (I14j|15[) . with three different sets 
of the particle numbers, N, and the lattice sizes, L: 
(N, L) = (4,4), (4,5), (5,5), which result in the Hilbert 
spaces of size, N$ — 35, 70, 126, respectively. The 
other parameters of the corresponding Hamiltonians are 
the same as in the previous section, i.e. J/U = 7/3, 
k B T A = 5{U + J), and k B T B = U + J. 

The numerical results demonstrate the presence of a 
smooth crossover between the arithmetic-mean and ther- 
mal regimes of equilibration. The regimes reveal them- 
selves by the plateau-like minima of X>p T ' can . The in- 
crease of system size Ms sharpens the crossover, note 
the log-scale in Fig. @] In addition, there are three no- 
table features: (i) the crossover region shifts to the re- 
gion of smaller coupling strengths, thus squeezing the 
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FIG. 5: The same distances as in Fig. [4J but as func- 
tions of the coupling strength rescaled by the mean 
level spacing of the corresponding system. After the 
scaling, in contrast to Fig. f4J the systems of different size 
share the same crossover region (open circle). Note, that due 
to the logarithmic scale of :r-axis, the region of the thermal 
equilibration is much broader than that of the arithmetic- 
mean equilibration. 



region of arithmetic-mean equilibration with increasing 
Ms; (ii) the region of thermal equilibration extends 
and the distance 2?™ n decreases with increasing Ms; 
(iii) the finite region of thermal equilibration is followed 
by the regime where the the weak coupling condition, 
Eq. (6) of main manuscript, is violated, causing the quasi- 
equilibrium energy level populations to deviate signifi- 
cantly from the thermal populations. 

The density of states for the bosonic model, Eq. (ITU) , 
with fixed Hamiltonian parameters and the mean occu- 
pation number (n) — N/L ~ 1, possesses essentially 
identical shapes. Thus we expect that the mean en- 
ergy level spacing, s — (en s — e\)/{Ms — 1), sets the 
proper energy scale for the crossover region between the 
arithmetic-mean and thermal equilibrations for different 
system sizes, Ms- 

The results of the scaling depicted in Figure [5] sup- 
port this hypothesis. The inset of Fig. [4] shows that the 
mean energy level spacing scales like No , with the expo- 
nent £ ~ —0.8. This also provides us with the scaling 
of the arithmetic-mean region with increasing of the sys- 
tem size. The scaling exponent seems to depend on the 
mean lattice occupation number (n) and the ratio of the 
Hamiltonian parameters, J/U, only. 



DIFFERENT SYSTEMS IN CONTACT 

Here we present the example of the equilibration pro- 
cess between two quantum systems with different spectra 
and initial temperatures, but the same size of the Hilbert 
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FIG. 6: Relaxation pathways for the bosonic model 
with two different systems in contact. Both systems are 
initially prepared in Gibbs states at different temperatures, 
(a) The evolution of the mean system energies and (b) the 
corresponding temperatures of the 'hot' system A (with N = 
4 bosons on the lattice with L = 6 sites) and the 'cold' system 
B (with N = L — 5) are shown by the red and blue lines, 
respectively. The time is plotted in units of inverse mean level 
spacing s of the system B. (c) The energy level populations 
of both the systems are displayed at different moments of 
time, marked by the corresponding symbols in (a, b). The 
solid line corresponds to the Gibbs energy level populations 
at the temperatures evaluated from the temporal values of 
mean system energies. The Hamiltonian parameters for both 
the system and the initial system temperatures are the same 
as in the case of the two identical systems in main manuscript, 
i.e. T A = 5(17 + J) = 94.91 s and T B = U + J= 18.98 s, while 
the coupling constant only is five times larger, XU lnt = 9.5 s. 



space, A/a = Afs, see Fig. [5] The lack of the permuta- 
tion symmetry, A «-> B, leads to the absence of two-fold 
degeneracies in the composite system at A = 0, meaning 
that an infinitesimal coupling strength A > does not 
modify significantly the spectrum of composite system 
any more, in contrast to the previous case of two iden- 
tical systems in contact. Therefore different subsystems 
do not exhibit the regime of arithmetic-mean equilibra- 
tion. At very weak coupling, which would induce the 
arithmetic-mean equilibration for two identical systems, 
the 'thermal equilibration' between different systems oc- 
curs locally, i.e. within a number of independent clusters 
of energy levels only. In order to render the thermal equi- 
libration process between the subsystems of rather small 
sizes, Ma = Mb = 126, we had to choose a relatively 
large value of the coupling constant, which is close to the 
upper bound of the weak coupling limit, Eq. (6) of main 
manuscript. 
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